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Abstract 

Over the past years, various representation systems which sparsely approximate functions 
governed by anisotropic features such as edges in images have been proposed. We exemplarily 
mention the systems of contourlets, curvelets, and shearlets. Alongside the theoretical devel- 
opment of these systems, algorithmic realizations of the associated transforms were provided. 
However, one of the most common shortcomings of these frameworks is the lack of providing 
a unified treatment of the continuum and digital world, i.e., allowing a digital theory to be a 
natural digitization of the continuum theory. In fact, shearlet systems are the only systems so 
far which satisfy this property, yet still deliver optimally sparse approximations of cartoon-like 
images. In this chapter, we provide an introduction to digital shearlet theory with a particular 
focus on a unified treatment of the continuum and digital realm. In our survey we will present 
the implementations of two shearlet transforms, one based on band-Umited shearlets and the 
other based on compactly supported shearlets. We will moreover discuss various quantitative 
measures, which allow an objective comparison with other directional transforms and an objec- 
tive tuning of parameters. The codes for both presented transforms as well as the framework 
for quantifying performance are provided in the Matlab toolbox jShearLab 

1 Introduction 

One key property of wavelets, which enabled their success as a universal methodology for signal 
processing, is the unified treatment of the continuum and digital world. In fact, the wavelet 
transform can be implemented by a natural digitization of the continuum theory, thus providing 
a theoretical foundation for the digital transform. Lately, it was observed that wavelets are 
however suboptimal when sparse approximations of 2D functions are sought. The reason is 
that these functions are typically governed by anisotropic features such as edges in images or 
evolving shock fronts in solutions of transport equations. However, Besov models - which 
wavelets optimally encode - are clearly deficient to capture these features. Within the model of 
cartoon-like images, introduced by Donoho in lfT3l in 1999, the suboptimal behavior of wavelets 
for such 2D functions was made mathematically precise; see also Chapter [4 |. 

Among the various directional representation systems which have since then been proposed 
such as contourlets lfT2l . curvelets ||9l, and shearlets, the shearlet system is in fact the only one 
which delivers optimally sparse approximations of cartoon-like images and still also allows for 
a unified treatment of the continuum and digital world. One main reason in comparison to the 
other two mentioned systems is the fact that shearlets are affine systems, thereby enabling an 
extensive theoretical framework, but parameterize directions by slope (in contrast to angles) 
which greatly supports treating the digital setting. As a thought experiment just note that a 
shear matrix leaves the digital grid 1? invariant, which is in general not true for rotation. 

This raises the following questions, which we will answer in this chapter: 

(PI) What are the main desiderata for a digital shearlet theory? 
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(P2) Which approaches do exist to derive a natural digitization of the continuum shearlet the- 
ory? 

(P3) How can we measure the accuracy to which the desiderata from (PI) are matched? 

(P4) Can we even introduce a framework within which different directional transforms can be 
objectively compared? 

Before delving into a detailed discussion, let us first contemplate about these questions on 
a more intuitive level. 

1.1 A Unified Framework for the Continuum and Digital World 

Several desiderata come to one's mind, which guarantee a unified framework for both the con- 
tinuum and digital world, and provide an answer to (PI). The following are the choices of 
desiderata which were considered in EOl 1141 : 

• Parseval Frame Property. The transform shall ideally have the tight frame property, which 
enables taking the adjoint as inverse transform. This property can be broken into the fol- 
lowing two parts, which most, but not all, transforms admit: 

o Algebraic Exactness. The transform should be based on a theory for digital data in 
the sense that the analyzing functions should be an exact digitization of the contin- 
uum domain analyzing elements. 

o Isometry of Pseudo-Polar Fourier Transform. If the image is first mapped into a 
different domain - here the pseudo-polar domain -, then this map should be an 
isometry. 

• Space-Frequency-Localization. The analyzing elements of the associated transform should 
ideally be highly localized in space and frequency - to the extent to which uncertainty 
principles allow this. 

• Shear Invariance. Shearing naturally occurs in digital imaging, and it can - in contrast to 
rotation - be precisely realized in the digital domain. Thus the transform should be shear 
invariant, i.e., a shearing of the input image should be mirrored in a simple shift of the 
transform coefficients. 

• Speed. The transform should admit an algorithm of order 0{N'^\ogN) flops, where N'^ is 
the number of digital points of the input image. 

• Geometric Exactness. The transform should preserve geometric properties parallel to 
those of the continuum theory, for example, edges should be mapped to edges in transform 
domain. 

• Stability. The transform should be resilient against impacts such as (hard) thresholding. 

1.2 Band-Limited Versus Compactly Supported Shearlet Transforms 

In general, two different types of shearlet systems are utilized today: Band-limited shearlet 
systems and compactly supported shearlet systems (see also Chapters HI and JU). Regarding 
those from an algorithmic viewpoint, both have their particular advantages and disadvantages: 

Algorithmic realizations of the band-limited shearlet transform have on the one hand typ- 
ically a higher computational complexity due to the fact that the windowing takes place in 
frequency domain. However, on the other hand, they do allow a high localization in frequency 
domain which is important, for instance, for handling seismic data. Even more, band-limited 
shearlets do admit a precise digitization of the continuum theory. 

In contrast to this, algorithmic realizations of the compactly supported shearlet transform 
are much faster and have the advantage of achieving a high accuracy in spatial domain. But for 
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a precise digitization one has to lower one's sights sHghtly. A more comprehensive answer to 
(P2) will be provided in the sequel of this chapter, where we will present the digital transform 
based on band-limited shearlets introduced in 1,20 J and the digital transform based on compactly 
supported shearlets from 1221 . 

1.3 Related Work 

Since the introduction of directional representation systems by many pioneer researchers (fSl 
|9, 10, 11 , 12|), various numerical implementations of their directional representation systems 
have been proposed. Let us next briefly survey the main features of the two closest to shearlets, 
which are the contourlet and curvelet algorithms. 

• Curvelets fT\. The discrete curvelet transform is implemented in the software package 
CurveLab, which comprises two different approaches. One is based on unequispaced 
FFTs, which are used to interpolate the function in the frequency domain on different 
tiles with respect to different orientations of curvelets. The other is based on frequency 
wrapping, which wraps each subband indexed by scale and angle into a fixed rectangle 
around the origin. Both approaches can be realized efficiently in 0{N^\ogN) flops with 
A' being the image size. The disadvantage of this approach is the lack of an associated 
continuum domain theory. 

• Contourlets |12|. The implementation of contourlets is based on a directional filter 
bank, which produces a directional frequency partitioning similar to the one generated 
by curvelets. The main advantage of this approach is that it allows a tree-structured filter 
bank implementation, in which aliasing due to subsampling is allowed to exist. Conse- 
quently, one can achieve great efficiency in terms of redundancy and good spatial local- 
ization. A drawback of this approach is that various artifacts are introduced and that an 
associated continuum domain theory is missing. 

Summarizing, all the above implementations of directional representation systems have their 
own advantages and disadvantages; one of the most common shortcomings is the lack of pro- 
viding a unified treatment of the continuum and digital world. 

Besides the shearlet implementations we will present in this chapter, we would like to refer 
to Chapter |2 | for a discussion of the algorithm in 1 16] based on the Laplacian pyramid scheme 
and directional filtering. It should be though noted that this implementation is not focussed on a 
natural digitization of the continuum theory, which is a crucial aspect of the work presented in 
the sequel. We further would like to draw the reader's attention to Chapter |3 1 which is based on 
II2TI aiming at introducing a shearlet MRA from a subdivision perspective. Finally, we should 
mention that a different approach to a shearlet MRA was recently undertaken in ifTTll . 

1.4 Framework for Quantifying Performance 

A major problem with many computation-based results in applied mathematics is the non- 
availability of an accompanying code, and the lack of a fair and objective comparison with 
other approaches. The first problem can be overcome by following the philosophy of 're- 
producible research' |15| and making the code publicly available with sufficient documenta- 
tion. In this spirit, the shearlet transforms presented in this chapter are all downloadable from 
|http : / /www . shearlab . org| One approach to overcome the second obstacle is the provi- 
sion of a carefully selected set of prescribed performance measures aiming to prohibit a biased 
comparison on isolated tasks such as denoising and compression of specific standard images 
like 'Lena', 'Barbara', etc. It seems far better from an intellectual viewpoint to carefully de- 
compose performance according to a more insightful array of tests, each one motivated by a 
particular well-understood property we are trying to obtain. In this chapter we will present 
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such a framework for quantifying performance specifically of implementations of directional 
transforms, which was originally introduced in f20','l4l. We would like to emphasize that such 
a framework does not only provide the possibility of a fair and thorough comparison, but also 
enables the tuning of the parameters of an algorithm in a rational way, thereby providing an 
answer to both (P3) and (P4). 



1.5 Shear Lab 



Following the philosophy of the previously detailed thoughts, |ShearLalq | was introduced by 
Donoho, Shahram, and the authors. This software package contains 



• An algorithm based on band-limited shearlets introduced in 

• An algorithm based on compactly supported separable shearlets introduced in 

• An algorithm based on compactly supported nonseparable shearlets introduced in 1*231. 

• A comprehensive framework for quantifying performance of directional representations 
in general. 

This chapter is also devoted to provide an introduction to and discuss the mathematical founda- 
tion of these components. 

1.6 Outline 

In Section|2] we introduce and analyze the fast digital shearlet transform FDST, which is based 
on band-limited shearlets. Section [3] is then devoted to the presentation and discussion of 
the digital separable shearlet transform DSST and the digital non-separable shearlet transform 
DNST. The framework of performance measures for parabolic scaling based transforms is pro- 
vided in Section]?] In the same section, we further discuss these measures for the special cases 
of the three previously introduced transforms. 

2 Digital Shearlet Transform using Band-Limited Shearlets 

The first algorithmic realization of a digital shearlet transform we will present, coined Fast 
Digital Shearlet Transform (FDST), is based on band-limited shearlets. Let us start by defining 
the class of shearlet systems we are interested in. Referring to Chapter lUl, we will consider 
the cone-adapted discrete shearlet system SH{(j),\(f, v/;A, A,A) — 4>(0;A) U'P(v/;A)U^(v/;A) 
with A = and 

A = A = {{j,k,m) : j>0,\k\ <2^,m(^l}}. 

We wish to emphasize that this choice relates to a scaling by 4^ yielding an integer valued 
parabolic scaling matrix, which is better adapted to the digital setting than a scaling by 2'. We 
further let be a classical shearlet (i// Ukewise with iff {^1,^2) = Wi^2,^i)), i-e., 

V/(<^) = V>(<^i,^2) = V>i(^i)r2(|), (2.1) 

where i/i e L^(IR) is a wavelet with i/i e C°°(M) and supp i/i C [— 4, U [5, 4], and 1/^2 € 
L^(]R) a 'bump' function satisfying \j/2 G C°°(K) and supp \j/2 C [—1,1]. We remark that the 
chosen support deviates slightly from the choice in the introduction, which is however just a 
minor adaption again to prepare for the digitization. Further, recall the definition of the cones 
■^11 - "^22 from Chapter Q. 



'ShearLab (Version 1.1) is available from http : / /www . shearlab . org 
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The digitization of the associated discrete shearlet transform will be performed in the fre- 
quency domain. Focussing, on the cone '€21, say, the discrete shearlet transform is of the form 



= (/,V>r,) = {f.2-n-^[SlA,-y)e^<^^->'^^--)), (2.2) 

where r\ — {j,k,m,i) indexes scale j, orientation k, position m, and cone i. Considering this 
shearlet transform for continuum domain data (taking all cones into account) implicitly induces 
a trapezoidal tiling of frequency space which is evidently not cartesian. A digital grid perfectly 
adapted to this situation is the so-called 'pseudo-polar grid', which we will introduce and dis- 
cuss subsequently in detail. Let us for now mention that this viewpoint enables representation 
of the discrete shearlet transform as a cascade of three steps: 

1) Classical Fourier transformation and change of variables to pseudo-polar coordinates. 

2) Weighting by a radial 'density compensation' factor. 

3) Decomposition into rectangular tiles and inverse Fourier transform of each tiles. 
Before discussing these steps in detail, let us give an overview of how these steps will be 



faithfully digitized. First, it will be shown in Subsection 2.1 that the two operations in Step 1) 
can be combined to the so-called pseudo-polar Fourier transform. An oversampling in radial 
direction of the pseudo-polar grid, on which the pseudo-polar Fourier transform is computed, 
will then enable the design of 'density-compensation-style' weights on those grid points leading 



to Steps 1) & 2) being an isometry. This will be discussed in Subsection 2.2 Subsection 2.3 
is then concerned with the digitization of the discrete shearlets to subband windows. Notice 
that a digital analog of ( |2.2| l moreover requires an additional 2D-iFFT. Thus, concluding the 
digitization of the discrete shearlet transform will cascade the following steps, which is the 
exact analogy of the continuum domain shearlet transform ( |2.2| l: 

(51) PPFT: Pseudo-polar Fourier transform with oversampling factor in the radial direction. 

(52) Weighting: Multiplication by 'density-compensation-style' weights. 

(53) Windowing: Decomposing the pseudo-polar grid into rectangular subband windows with 
additional 2D-iFFT. 

With a careful choice of the weights and subband windows, this transform is an isometry. Then 
the inverse transform can be computed by merely taking the adjoint in each step. A final dis- 



cussion on the FDST will be presented in Subsection 2.4 



2.1 Pseudo-Polar Fourier Transform 

We start by discussing Step (SI). 



2.1.1 Pseudo-Polar Grids with Oversampling 

In pSl, a fast pseudo-polar Fourier transform (PPFT) which evaluates the discrete Fourier trans- 
form at points on a trapezoidal grid in frequency space, the so-called pseudo-polar grid, was 
akeady developed. However, the direct use of the PPFT is problematic, since it is - as de- 
fined in |5 | - not an isometry. The main obstacle is the highly nonuniform arrangement of 
the points on the pseudo-polar grid. This intuitively suggests to downweight points in regions 
of very high density by using weights which correspond roughly to the density compensation 
weights underlying the continuous change of variables. This will be enabled by a sufficient 
radial oversampling of the pseudo-polar grid. 
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This new pseudo-polar grid, which we will denote in the sequel by Q.r to indicate the over 
sampling rate R, is defined by 



(2.3) 



where 



tl r'N'r)- 2-^-2' 2 -"-2 J' 



2n le 
' R ' N 



2 * 2 ' 
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RN 



}■ 



(2.4) 
(2.5) 



This grid is illustrated in Fig. [T] We remark that the pseudo-polar grid introduced in 151 coin- 
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Figure 1 : The pseudo-polar grid Q.r = Q.\ U iorN = 4 and 7? = 4. 

cides with Q.r for the particular choice R = 2. It should be emphasized that Q.g = U £2| is 
not a disjoint partitioning, nor is the mapping (n,£) H> • ^, ^) or — ^ • ^) injective. 
In fact, the center 

■^-{(0,0)} (2.6) 
appears + 1 times in Q.^ as well as £2^, and the points on the seam lines 

-^R = <«<f ,«^0}, 

'^R = {(|,-|):-f <«<f ,«^0}. 

appear in both ilj^ and Q.^. 

Definition 2.1. Lef A^, R be positive integer, and let D,r be the pseudo-polar grid given by \23\ . 
For an N X N image I := {/(m, v) ■ < «, v < ^ — 1}, the pseudo-polar Fourier transform 
(PFFT) I of I evaluated on Q,r is then defined to be 

N/2-1 j^. 

7(0)1,0)2)= /(«,v).-^(""'+'"^\ (0)1,0)2) Gn«, 

ii,v=-N/2 

where nn)>N is an integer. 

We wish to mention that mg > is typically set to be niQ = ^{RN + 1) for computational 
reasons (see also [Sj), but we for now allow more generality. 

2.1.2 Fast PPFT 

It was shown in jS), that the PPFT can be realized in 0{N'^\ogN) flops with N x N being the 
size of the input image. We will now discuss how the extended pseudo-polar Fourier transform 
as defined in Definition |2.1| can be computed with similar complexity. 
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For this, let / be an image of size N x N. Also, mo is set - but not restricted - to be 
mo = ^{RN + 1); we will elaborate on this choice at the end of this subsection. We now focus 
on Q.^,, and mention that the PPFT on the other c one can be computed similarly. Rewriting 
the pseudo-polar Fourier transform from Definition 2 . 1 for (0)1,0)2) = (— ^ ■ j^, ^) € , we 
obtain 

N/2-1 . . N/2-1 N/2-1 , . , ,, , 

7(0)1,0)2) = 2^ I{u,v)e "'0^ 1^ 1^ I{u,v)e -"o ^ 

u.y=-N/2 u=-N/2v=-N/2 
N/2-1 ( N/2-\ \ , 

= 2^ 2- Ku,v)e nN+i e (av+i)jv_ (2.7) 

u=-N/2 \v=-N/2 ) 



This rewritten form, i.e., p.7| l, suggests that the pseudo-polar Fourier transform 7 of / on £2^ 
can be obtained by performing the ID FFT on the extension of 7 along direction v and then 
applying a fractional Fourier transform (frFT) along direction u. To be more specific, we require 
the following operations: 

Fractional Fourier Transform. For c G C^^', the (unaliased) discrete fractional Fourier 
transform by a ^ C is defined to be 



N/2 

iFff^,c){k):= £ c(y).-2--->^-«, fc^-f,...,f. 



= -N/2 

It was shown in {SI, that the fractional Fourier transform F^^yC can be computed using 0{N\ogN) 
operations. For the special case of a = 1/(A^ + 1), the fractional Fourier transform becomes 
the (unaliased) ID discrete Fourier Transform (ID FFT), which in the sequel will be denoted 
by 7^1. Similarly, the 2D discrete Fourier Transform (2D FFT) will be denoted by 7^2, and the 
inverse of the 72 by 7^2^' (2D iFFT). 

Padding Operator. For even, m > N an odd integer, and c G C'^, the padding operator 
Em.n gives a symmetrically zero padding version of c in the sense that 



Using these operators, \2.1\ can be computed by 

N/2-1 



I{cOi,0>2) = FioERN+i_NoI{u,n)e 

u=-N/2 
N/2 



-2ttuiP- — 



= L EN+x.NoF,oERN+x.NoI{u,n)e 2™^' 

u=-N/2 



= iF;}'uii;n)m, 

E««+L«o7GC(«^+^) 

and ID frFT require only 0{NlogN) operations for a vector of size A^, the total complexity of 



where7 = £A,+i,A,oFio£KAr+i^A,o7GC(^"+^'^('^+i) and a„ = - p^^^^ . Since the ID FFT 



this algorithm for computing the pseudo-polar Fourier transform from Definition 2.1 is indeed 
0{N^\ogN) for an image of size N xN. 

We would like to also remark that for a different choice of constant mo, one can compute the 
pseudo-polar Fourier transform also with complexity 0{N-^logN) for an image of size N xN. 
This however requires application of the fractional Fourier transform in both directions u and v 
of the image, which results in a larger constant for the computational cost; see also [6 |. 
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2.2 Density-Compensation Weights 



Next we tackle Step (S2), which is more deUcate than it might seem, since the weights will not 
be derivable from simple density compensation arguments. 



2.2.1 A Plancherel Theorem for the PPFT 

For this, we now aim to choose weights w : CIr M+ so that the extended PPFT from Definition 



2.1 becomes an isometry, i.e.. 



N/2-l 

^ |/(«,V)P= E Vv(a)l,C02)- 1/(0)1, 0)2) |2. 



(2.8) 



Observing the symmetry of the pseudo-polar grid, it seems natural to select weight functions w 
which have full axis symmetry properties, i.e., for all {(Oi , (02) G ^r, we require 



w{C0i,ah) = w{C02,C0i), w{(0i,0}2) = w{-C0i,(J>2), ^(0)1,0)2) = w{C0i,-(J>2). 



(2.9) 



Then the following 'Plancherel theorem' for the pseudo-polar Fourier transform on £2^ - similar 
to the one for the Fourier transform on the cartesian grid - can be proved. 

Theorem 2.1 (I120j). Let N be even, and let w : £2^ — > be a weight function satisfying ( |2.9| ). 
Then ( |2.8| l holds, if and only if, the weight function w satisfies 

5{u,v) = w(0,0) 

RN/2 

+ 4- E 1: w(|,|.-^).cos(2;r«.|^).cos(2;rv.^.^) 

e=0,N/2 n=l 
N/2-1RN/2 

+ 8-1 1: w(|,|.-^).cos(27r«.^).cos(2;rv.^.|) (2.10) 

e=i n=l 

forall-N+l <u,v<N-l. 

Proof. We start by computing the right hand side of ( |2.8| ): 

E Vv(t(Ji,C(}2)- |/(«l, 0)2)1^ 
{(Ul,(U2)Gi2fi 



E w(a)i,a)2)' 

E w{C0u0>2)' 



E w(a)i,a)2)- E |/(m,v)| 

{mi,o>2)eQ.R u,v=-N/2 



N/2-1 

E /(«,v).-^("'"'+^"^' 

u,v=-N/2 

N/2-1 N/2-1 , , ■ 

E E /(«,v)7(«vo.-^«"-""*"+(''-^'^'^' 

_u,v=-N/2u'y=-N/2 
N/2-1 

2 



N/2-1 



£ /(«,v)/(«',v'). 



u,v,u',v'=-N/2 
(m,v)^(«',v') 



52 w(a)i,a)2)-e 



-||((«-«')'«l+(v-v')'i>2) 



Choosing / = c„,,vi5(m — mi;V — vi) + Cu^^y^5{u — U2,v — V2) for all —N/l < mi, vi, M2, v'2 < 
N/2 — 1 and for all c„, vi jCu^ y, G C, we can conclude that p.8| l holds if and only if 

£ w(tOi,«2)-e"^*"'"+""^'=5(M,v), -A?+l <M,v<A?-l. 

((Bl,(02)Gn« 
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By the symmetry of the weights ( |2.9| l, this is equivalent to 

^ w{(OuO>2)-[cos{^u(Oi)cos{^vco2)]^5{u,v) (2.11) 

for all — 1 <M,v<A^— 1. From this, we can deduce that ( |2.1 1[ ) is equivalent to ( |2.10[ ), 
which proves the theorem. □ □ 

Notice that ( |2.10| l is a lineai" system with RN^/4 + RN/2 + 1 unknowns and {2N - 1)^ 
equations, wherefore, in general, one needs the oversampling factor R to be at least 16 to enforce 
solvability. 



2.2.2 Relaxed Form of Weight Functions 



The computation of the weights satisfying Theorem 2.1 by solving the full linear system of 
equations ( |2.10| i is much too complex. Hence, we relax the requirement for exact isometric 
weighting, and represent the weights in terms of undercomplete basis functions on the pseudo- 
polar grid. 

More precisely, we first choose a set of basis functions wi , . . . , w„q : D.R — > M+ such that 

"0 

Y^Wjicou (02)^0 for all (0)1,0)2) e 
We then represent weight functions w : Q.r by 

no 

w:=^CjW^-, (2.12) 

;=i 

with ci,....c„Q being nonnegative constants. This approach now enables solving ( |2.10| i for 
the constants ci,...,c„g using the least squares method, thereby reducing the computational 
complexity significantly. The 'full' weight function w is then given by ( |2.12[ ). 

We next present two different choices of weights which were derived by this relaxed ap- 
proach. Notice that (o)i , 0)2) and {n,£) will be used interchangeably. 

Choice 1. The set of basis functions wi , . . . , W5 is defined as follows: 
Center. 

= 1(0.0). 

Boundary: 

^2 — l{((i)|,(i)2):|n|=Afff/2, 01=0)2} ™d W3 ~ 1 {(o, ,(i>2):|n|=Arfi/2, (Oi /fih} i 

Seam lines: 

W4 = \n\ ■ l{(ffli,o)2):l<|n|<A'«/2,(Bi=(B2}. 

Interior. 

^5 = \n\ ■ l{((ai,(B2):l<|n|<Af«/2,e),/(i)2}- 

Choice 2. The set of basis functions wi , . . . , /2+2 is defined as follows: 
Center. 

= 1(0^0). 

Radial Lines: 

^^+2 = l{(o,,O2):l<|„|<™/2,0>, = 45<i„}' ^ = 0, 1, . . . ,N/2. 

The associated weight functions are displayed in Fig. [2] In general, suitable weight func- 
tions usually obey the pattern of linearly increasing values along the radial direction. Thus, this 
is a natural requirement for the basis functions. 
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Choice 1 



Choice 2 



Figure 2: Weight functions on the pseudo-polar grid for N = 256 and 7? = 8. 



2.2.3 Computation of the Weighting 



For the FDST - as also in the implementation in ShearLab'- the coefficients in the expansion 
p.l2[ ) will be computed off-line, and then hardwired in the code. This enables the weighting of 
a function on the pseudo-polar grid to simply be a point-wise multiplication in each sampling 
point. That is, letting 7 ; = / : Q.R C be the pseudo-polar Fourier transform ofanNxN image 
/ and w : Q.r be any suitable weight function on £1/;, the values 

A,(«i,C02) (02)- v^vv(0Ji7ft^ for all (couCOi) e 

need to be computed. 

Let us comment on why the square root of the weight is utilized. If the weights w satisfy 
the condition in Theorem 2.1 we obtain P*wP = Id, which can be written in a symmetric form 
as follows: {y^P)* y/wP = Id. This form shows that the operator y^P can be inverted by 
taking the adjoint {y/wP)*. In other words, each image can be reconstructed from its weighted 
pseudo-polar Fourier transform by applying the adjoint of the weighted pseudo-polar Fourier 
transform. This issue will be discussed in further detail in Subsection l2.4.2l 



2.3 Digital Shearlets on Pseudo-Polar Grid 

We next aim at deriving a faithful digitization of the shearlet transform associated with a band- 
limited cone-adapted discrete shearlet system to the pseudo-polar grid. This would settle Step 
(S3). 



2.3.1 Preparation for Faithful Digitization 

For this, let us recall the definition of the discrete shearlet transform associated with ( |2.2| l; taking 
the particular form ( |2.1| ) of the shearlet Xj/ G L?{M?) into account. Restricting our attention to 
the cone '^21^ we obtain 

/ ^ (/,2-4v>(5[A4-r)Z^,,e2-(V.^*-'->^ 

= (/, 2-^5 (4-^'<^i ) V>2 {k + 2' I )xy,,, e2-(^4-.^*'"-> ) , 

for scale j, orientation k, position m, and cone i. 

To approach a faithful digitization, we first have to partition Q.r according to the partitioning 
of the plane into '^^n, ^€\2, '^2\, and ^ri, as well as a centered rectangle The center as 
defined in ( |2.6| l will play the role of S^,. Thus it remains to partition the set Q.r beyond the 
already defined partitioning into £1^ and n| (cf. ( |2.4| i and ( |2.5| l) by setting 

Q^, = Q}}y^'gy^Q}K and n| = n^'u^unF, 
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where 





= {(-f 


2<; 


2n\ 
R ) 


N 
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<e< 


N 
2 


,l<n<f }, 




q12 


= {(-f 


2^ 
N 


2n\ 
R ) 


N 
2 


<i< 


N 
2 


<n<- 


1} 




= 


In 
R 


^) 

ni 


N 
2 


<i< 


N 
2 


,l<n<f} 




iZj; 


= 


In 
R 


^) 


N 
2 


<(■< 


N 
2 


<n<- 


1} 



When restricting to the cone Q.^ , say, the exact digitization of the coefficients of the discrete 
shearlet system is 

(B:=((Bi,(B2)en|i 

£ 7(fi)i,a)2)2-^iw(4->a),)y(A: + 2^'^)e-2'^'(Vi^^'«''») 

{coi.,m2)enl^ 



RN 



£ £ /(©i,fi>2)2-^'2W(4-^'^)y(fe-2^'+i|)e-2'^'('"'^^Vi»), (2.13) 



n=l / 



where V and W as well as the ranges of j, k, and m are to be carefully chosen. 

Our main objective will be to achieve a digital shearlet transform, which is an isometry. 
This - as in the continuum domain situation - is equivalent to requiring the associated shearlet 
system to form a tight frame for functions / : CIr C. For the convenience of the reader let us 
recall the notion of a Parseval frame in this particular situation. A sequence {(pi)xeA - ^ being 
some indexing set - is a tight frame for all functions J : Q.R ^ C, ii 



E 

XeA 



£ 7(C0i,a)2)(PA(t»l,t»2) 

((»l,(02)GnR 



E IA«1,»2)P. 

{u>i.u>2)eSlR 



In the sequel we will define digital shearlets on Q.^^ and extend the definition to the other 
cones by symmetry. 



2.3.2 Subband Windows on the Pseudo-Polar Grid 

We start by defining the scaling function, which will depend on two functions Vq and Wo, and 
the generating digital shearlet, which will depend on again two functions V and W. Wq and W 
will be chosen to be Fourier transforms of wavelets, and Vq and V will be chosen to be 'bump' 
functions, paralleling the construction of classical shearlets. 

First, let Wq be the Fourier transform of the Meyer scaling function such that 

supp Wo C [-1, 1] and Wb(±l) = 0, (2.14) 

and let Vq be a 'bump' function satisfying 

supp VqC [-3/2,3/2] with yo(§) = 1 for |^| < 1,^ e K. 

Then we define the scaling function (j) for the digital shearlet system to be 

^{^u^2)=Woi4-J^^i)Vo{4-J^^2), (^1,^2) eM'. 

For now, we define it in continuum domain, and will later restrict this function to the pseudo- 
polar grid. 
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Let next W be the Fourier transform of the Meyer wavelet function satisfying the support 
constraints 

supplVC [-4,1/4] U [1/4,4] and lV(±l/4) = 1V(±4) = 0, (2.15) 
as well as, choosing the lowest scale ji to be ji := — \log^{R/2)~\ , 

riog4A'i 

\Wq{4-J''^)\^+ £ |W(4-^'^)|2 = l forall 1^1 <A?, ^ eM. (2.16) 
;=./L 

We further choose V to be a 'bump' function satisfying 

suppyc[-l,l] and y(±l)=0, (2.17) 

as well as 

|y(^-l)|2 + |y(^)|2 + |y(^ + l)|2 = i forall|^|<l,^eM. (2.18) 
Then the generating shearlet Xj/ for the digital shearlet system on il^ is defined as 

v/(<^i,&)=w(<^,)y(|), {^u^2)eR\ (2.19) 

Notice that p.l8| l implies 

2' 

£ |y(2>^ -s)P = 1 forall 1^1 < 1, <^ eM;; >0, (2.20) 

which will become important for the analysis of frame properties. For the particular choice of 
Vo, Wo, V, and W injShearLab) we refer to Subsection |2.3.7| 

2.3.3 Range of Parameters 

We from now on assume that R and are both positive even integers and that — 2"" for some 
integer no G N. This poses no restrictions, since both parameters can be enlarged to satisfy this 
condition. 

We start by analyzing the range of j. Recalling the definition of the shearlet \j/ in ( |2.19[ l 
and the support properties of W and V in p.l5| l and ( |2.17| l, respectively, we observe that the 
digitized shearlet 

2->iw(4--''f )y(/fc-2''+i|)e2'^'('"^^tV^'") (2.21) 
from ( |2.13[ ) has radial support 

n = 4'-'f+fi, f, =0,..., 4^-1 -if (2.22) 

on the cone Q.^ . To determine the appropriate range of j, we will analyze the precise support 
in radial direction. If j < — \log{R/2)~\ , then n < 1, which corresponds to only one point - the 
origin - and is dealt with by the scaling function. If j > [log4A^], we have n > Hence 
the value W{1 /4) = (cf. p.l5[ )) is placed on the boundary, and these scales can be omitted. 
Therefore, the range of the scaling parameter will be chosen to be 

j e {Jl, ■■■Jh}, where jt := - \\og{R/2)~\ and jn := [log4A^] . 

Next, we determine the appropriate range of k. Again recalling the definition of the shearlet 
yf in ( |2.19[ ), the digitized shearlet p.21| i has angular support 

e^2-j-^Nik-l)+t2, t2^0,...,2-jN (2.23) 
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on the cone £2|^ . To compute the range of k, we start by examining the case j >0. If k > 2^, 
we have £ > N/2. Hence the value 1) = (cf. ( |2.17[ )) is placed on the seam line, and these 
parameters can be omitted. By symmetry, we also obtain k > —iK Thus the shearing parameter 
will be chosen to be 

k&{-V,...,2^}. 

2.3.4 Support Size of Shearlets 

We next compute the support as well as the support size of scaled and sheared version of digital 
shearlets. This will be used for the normalization of digital shearlets. 

As before, we first analyze the radial support. By ( |2.22| i, the radial supports of the windows 
associated with scales Jl < j < Jh is 

„^4^-i|+fi, fj^o,...,4^-^-if , (2.24) 

and the radial support of the windows associated with the scale y'^ = —\log^{R/2)~\ and /// ~ 
[log4A^] are 

n = tu ri = l,...,4>-+lf, for; = ji, 

„=4^«-i|+fi, fi=o,...,f -4^«-'f, for; = ;^. 

Turning to the angular direction, by ( |2.23| l, the angular support of the windows at scale j 
associated with shears — 2-' < k <2-' is 

e^2-j-^N{k~l)+t2, t2 = 0,...,2-'N, (2.26) 

the angular support at scale j associated with the shear parameter k — —2^ is 

e = 2-j-^N{-2j - 1) +f2, t2=2-j^,.. .,2-jN, 

and for k — 2^ it is 

£ = 2-j-^N{2j -l)+t2, f2-0,...,2-^'f . (2.27) 

For the case j < 0, we simply let ^ = and £ — —N/2 + 12 with f? ~Q,.--,N. Also, for this 
lower frequency case, the window function W{4^^COi)V{k + 2^^) is slightly modified to be 

Wi4-jcoi)Vo{k + 2J^). 

These computations now allow us to determine the support size of the function W (4^^ a)i)y(fc - 
2-'^) in terms of pairs {n,£), which for scale j and shear k, is 



^l = J 4J-i.lf + i 
f-4'-if + l 



jL<j<jH, (2.28) 

J = Jh, 



and 




-2j <k<2j with > 0, 

k e {-V,V} with /■ > 0, (2.29) 

;<0. 



2.3.5 Digitization of the Exponential Term 

We next digitize the exponential term in ( |2.21| i, which can be rewritten as 



We observe two obstacles: 
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• The change of variables T := SlA^-jCO possible in ( |2.13[ ) can not be performed similarly 
in this situation due to the fact that the pseudo-polar grid is not invariant under the action 
of S^A^-j. This is however the first step in the continuum domain reasoning for tightness; 
see Chapter m. 

• The Fourier transform of a function defined on the pseudo-polar grid does not satisfy any 
Plancherel theorem. 

These problems require a slight adjustment of the exponential term, which will be the only 
adaption we allow us to make when digitizing. This will circumvent the two obstacles and 
enable us to construct a Parseval frame as well as derive a direct application of the inverse Fast 
Fourier transform in FDST. 

The adjustment will be made by using the mapping : M \ {0} M. defined by 9{x,y) = 
(x, j). This yields the modified exponential term 

-2^i(,«,(eo(s[)-')(4-^t,4-^'^^-2-^'^)) ^g-27ri(m,(4-^-^.-2^+l|)) 



(2.30) 



which can be rewritten as 



,-2;r,(«.(4-^ f,-2./+ljf)) ^ ^-2;r/(!^+(l-i)m2) -2;r,-(M,(4--'-^^-2^+l |)) 



with f 1 and f2 ranging over an appropriate set defined by ( |2.24| i, ( |2.25| l, and p.26| l-( |Z!27] i. Fig.|3] 
illustrates this adjustment. 



is] 



e 



Figure 3: Adjustment of the exponential term through the map 6 o[S 



Now, taking into account of the support size of each W{4 C0i)V {k + 2^ ^) as given in 



0)1 

( |2.28| ) and ( |2.29[ ), we obtain the following reformulation of ( |2.30| l: 

exp|-27i:i( m, ( -^-^T fl, — '--^^ j/J' ' 

This version shows that we might regard the exponential terms as characters of a suitable locally 
compact abelian group (see 1 18 1) with associated annihilator identified with the rectangle 

jl^T^,-^^) :n=0,...,if}^-l,r2=0,...,J-?,-l 

where if/ and ^f ,^ were defined in ( |2.28[ ) and ( |2.29| l, respectively. This viewpoint will be cru 




cial to guarantee that the digital shearlet system defined in Subsection 2.3.6 provides a Parseval 
frame on the pseudo-polar grid Q.r. In practice, ( |2.31| l also ensures that in Step (S3) on each 
windowed image on the pseudo-polar grid only a 2D-iFFT - in contrast to a fractional Fourier 
transform - needs to be performed, thereby reducing the computational complexity. 
For the low frequency square, we further require the set 

^ = {(n,r2):n = -l,...,l,r2--f,...,f}, 

which will be shown to be sufficient for guaranteeing that digital shearlet system forms a Par- 
seval frame. 
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2.3.6 Digital Shearlets 

We are now ready to define digital shearlets, which we define as functions on the pseudo-polar 
grid Q.j(. The spatial domain picture can thus be derived by the inverse pseudo-polar Fourier 
transform. 



Definition 2.2. Retaining the definitions and notations from Subsection \2.3\ for all (oi, Oi) G 
we define digital shearlets at scale j e {j^, . . . Jh}, shear k — [—2-', 2-'] H Z, and spatial 
position m € ^j,k by 



t21 



{0)1,(02) 



W {4-jcoi ) {k + 2j f )Xa2i (fi)i , 0)2) e 



27Ci{m.(4-J(Oi.2'^)) 



where — V for j >0 and = Vofor j < 0, and 

{couco2)^y^yjyl 

C{couco2) = { ti ■ («i,c%)e(^iu^|)\'^, 

on the remaining cones are defined accordingly by symmetry 
with equal indexing sets for scale j, shear k, and spatial location m. For Iq = 1,2, (01,0)2) G 
n]^', and no G we define the scaling function 



The shearlets (7, 
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.12 



^22 



j,k,m^ j,k^m^ j,k,m 



(p^ ((01,0)2) = 



C{ii>i.m2) 



H(Ol,0}2)Xa'o{cOh(02)e 



27r!-(no.(f.jvTl 



)) 



Then the digital shearlet system DSH is defined by 

DSH = : lo = l,2,«o e ^}U {(j|,,„ : ; G {Jl, • • ■,jH},ke {-2', 2'}, 

m€^^>,l = 11,12,21,22}. 

As desired, the digital shearlet system DSH, which we derived as a faithful digitization of 
the continuum domain band-hmited cone-adapted discrete shearlet system, forms a Parseval 
frame for J : Q.R C. 



Tlieorem 2.2 (120|). The digital shearlet system DSH defined in Definition 2.2 forms a Parseval 
frame for functions J : — > C. 



Proof. Letting J : £2^ C, we claim that 



'0,no lj,k,m 



(2.32) 



which proves the result. Here (7i,72)nR := L((Ui (i)2)GnR-^i(®i'^)-^2(«i,C^) for7i,72 : ^ 
C. 

We start by analyzing the first term on the RHS of ( |2.32| l. Let Iq G { 1 , 2} and Jc : ^r — > C be 
defined by Jc{(0i,(02) := C(cOi, 0)2) 7(0)1,0)2) for (o)i, 0)2) G £2^. Using the support conditions 
of 



7(0)1, 0)2)(p^° (0)1,0)2) 



"0 



7c(0)i, 0)2) • 0(0)1,0)2) 



. , N/2 

L L 7c(o)i, 0)2) -0(0)1, 0)2) -.-^-V' 



?';vtt)) 



"0 n=-ie=-N/2 
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The choice of ^ now allows us to use the Plancherel formula, see Subsection 2.3.5 Exploiting 
again support properties (see Subsection |2.3.5| l, we conclude that 

Ll(-/,1>^°)a,|'= E |C(«l,C02)-y(«i,C02)|2. 10(0)1, 0)2)^. 



Combining Iq = 1,2 and using ( |2.14| i, we proved 



EEl(-/,Onj'= E \Jicoi,a>2)\'-\Wo{co,)\\ 



(2.33) 



'0 «0 



Next we study the second term on the RHS in ( |2.32| l. By symmetry, it suffices to consider 
the case i =21. By the support conditions on W and V (see ( |2.15[ ) and ( |2.17[ )), 



E l(-^'^'UnJ' = E E E Jicouc^)cyjl„icouco2) 

hk.m j.kme.'if'jj, (cdi ,(U2)Ga|.' 



Lu^i E E JcicouO}2)-W{4-J(ai) 

j,k \-^jM ,neMj,t (fl),,(02)en|' 

J 4^'+l(«/2) 2--'-'A'(<r+l) 

Ltw^\ EE E ■/c(«i,o)2) 



Similarly as before, the choice of ^ does allow us to use the Plancherel formula, see Sub- 
section |233] Hence, 



Next we use ( |2.20[ ) to obtain 



E E ^c(«i,C02)-W(4-^a)i)-y^(^ + 2^-^ 



7/f 



E |7c(tOi,t02)|2 E |W^(4-^»i)l'- E l^''(^ + 2^' 

{ 61)1, 692 )Gn21 



J=.IL 



k=-V 



E \Jc{m,m)\^ E |w'(4-^«i)l'. 

{6l),,6l)2)Gn21 



]=]L 



Thus the second term on the RHS in ( |2.32| i equals 



E E l(-/>^kJn,i'= E E ■|w(4-^«i)i^ 

1 j-k.m (a\.a>2)eQ.R 



]=]L 



(2.34) 



Finally, our claim ( |2.32[ ) follows from combining ( |2.33| l, ( |2.34| i, and ( |2.16| ). □ □ 
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2.3.7 Digital Shearlet Windowing 



The final Step (S3) of the FDST then consists in decomposing the data on the points of the 
pseudo-polar grid given by the previously - in Steps (SI) and (S2) - computed weighted 
pseudo-polar image : Q.f; — > C into rectangular subband windows according to the digital 



shearlet system DSH defined in Definition 2.2 followed by a 2D-iFFT. More precisely, given 
Jw, the set of digital shearlet coefficients 



(■^»"^"o)a« for all lo, no 



and 



j,k.m ■ 



Jw. (J 



'i,k.m 



for all j,k,m, i 



is computed followed by application of the 2D-iFFT to each windowed image Jw(Pq and -/hC] ^ o 
restricted on the support of (p^ and ff^' ^ q, respectively. 

The definition of the digital shearlet system DSH in Definition 2.2 requires appropriate 



choices of the functions 0, Vq, V, Wq, and W, and the required conditions are stated throughout 
Subsection |2.3.2| We now discuss one particular choice, which is chosen injShearLabJ We 
start selecting the 'wavelets' Wq and W. In Subsection 2.3.2 these functions were defined to be 
Fourier transforms of the Meyer scaling function and wavelet function, respectively, i.e.. 







,[fv( 



\^\<l 

l<\^\<h 

otherwise, 



and 



sin[fv(l|^|-i)] 
W{^)^{ cos[fv(i|^|-i)] 




?<I^I<1, 
1<I^I<4, 
otherwise, 



where v > is a C*^ function or C°° function such that v(x) + v(l — x) = 1 for < x < 1. One 
possible choice for v is the function v{x) — x^(35 — S4x + 70x^ — 20x^), < x < 1, which 
then automatically fixes Wq and W. Since |Wo(^)P + |VK(i^)p = 1 for |^| < 1, the required 
condition ( |2.16[ ) is satisfied. The function v can be also used to design the 'bump' func- 
tion V as well, which needs to satisfy ( |2.18[ ). One possible choice for V is to define it by 
V{^) — \/v{l +1^) + v(l —^), — 1 < ^ < 1. Vb can then simply be chosen as Vb = 1. Let us 
finally mention that (j) is defined depending on Vo and Wq, wherefore fixing these two functions 
determines uniquely. 



2.4 Algorithmic Realization of the FDST 

We have previously discussed all main ingredients of the fast digital shearlet transform (FDST) 
- Fast PPFT, Weighting, and Digital Shearlet Windowing -, and will now summarize those 
findings. Depending on the application at hand, a fast inverse transform is required, which we 
will also detail in the sequel. In fact, we will present two possibilities: the Adjoint FDST and the 
Inverse FDST depending on whether the weighting allows to use the adjoint for reconstruction 
or whether an iterative procedure is required for higher accuracy. Fig.|4]provides an overview 
of the main steps of of the FDST and its inverse. For a more detailed description of FDST, 
Adjoint FDST, and Inverse FDST in form of pseudo-code, we refer to [20]. 

For the sake of brevity, we now let P, w, and W denote the Fast PPFT from Subsection |2.1.2| 



the weighting on the pseudo-polar grid described in Subsection 2.2.3 and windowing operator 



consisting of the application of the shearlet windows followed by 2D-iFFT to each array as 



detailed in Subsection 2.3.7 respectively 
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Figure 4: Flowcharts of the FDST (left) and its inverse (right). 



2.4.1 FDST 

We can summarize the steps of the algorithm FDST as follows: 



Step (SI): For a given image /, apply the Fast PPFT as described in Subsection 2.1.2 
obtain the function PI : Q.^ ^ C 

Step (S2): Apply the square root of an off-line computed weight function w : CIr 
PI as described in Subsection 



to 



to 



2.2.3 



yielding %/wPI : Q.r 



Step (S3): Apply the shearlet windows to the function wPI, followed by a 2D-iFFT to 
each array to ol 

ckm'i'^' '"pl- 



each array to obtain the shearlet coefficients W\/wPI, which we denote by c^", lo,no and 



2.4.2 Adjoint FDST 

Assuming that the weight function w used in Step (S2) satisfies the condition in Theorem |2.1| 
and using the Parseval frame property of the digital shearlet system (Theorem |2.2| i, we obtain 

{Wy/wP)*WV^P = P*Vw{W*W)V^P = P*wP = Id. 

Hence in this case, the FDST, which is abbreviated by Wy/wP can be inverted by applying the 
adjoint FDST, which cascades the following steps: 

• Step 1: For given shearlet coefficients C, i.e., c^q, io,no and c'-^^, j,k,m,i, compute 

the linear combination of the shearlet windows with coefficients cjjj,, lo,no and c^ j,,„, 
j,k,m,i. This gives the function W*C : Q.r — > C. 

• Step 2: Apply the square root of an off-line computed weight function w : Q.r — > C to 
W*C, yielding the function y^W*C : £2^ C. 

• Step 3: Apply the Fast Adjoint PPFT by running the Fast PPFT 'backwards' . For this, we 
just notice that the adjoint fractional Fourier transform of a vector c G C'^^' with respect 
to a constant a € C is given by F^^^c. Also, for m > A^, the adjoint padding operator 
£* ;v applied to a vector c e C" is given by {E*jjc){k) = c{k), k = ~N /2, . . . ,N /2 - 1. 
The Adjoint PPFT gives an image P*y/wW*C. 
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2.4.3 Inverse FDST 



Normally - as also with the relaxed form of weights debated in Subsection|2.2.2|- the weights 



will not satisfy the conditions of Theorem 2.1 precisely. A measure for whether application 



of the adjoint is still feasible will be discussed in Subsection 4.2 If higher accuracy of the 
reconstruction is required, one might use iterative methods, such as conjugate gradient methods. 
Since the digital shearlet system forms a Parseval frame, we always have 

W*W^/wP = V^P. 

Hence, iterative methods need to be 'only' applied to reconstruct an image / from knowledge 
of J := \fwPl, i.e., to solve the equation 

P*wPl = P*wJ 

for /. Since J might not be in the range of P, I is typically computed by solving the weighted 
least square problem min/ ||y^P/— -\/w-^I|2- Since the matrix corresponding to P*P is sym- 
metric positive definite, iterative methods such as the conjugate gradient method are applicable. 
The conjugate gradient method is then applied to the equation Ax = b with A ~ P*wP and 
b = P*wJ. Its performance can be measured by the condition number of the operator P*wP: 
cond{P*wP) = Xmax{P*wP)/Amin{P*wP), and it turns out that the weight function serves as a 
pre-conditioner. We remark that this measure is more closely studied in Subsection |4.2| 

To illustrate the behavior of the weights with respect to this measure, in Table[T]we compute 



cond{P*wP) for the weight functions arising from Choices 1 and 2 (cf. Subsection |2.2.2 1 with 
oversampling rate /? = 8. 

Table 1: Comparison of Choices 1 and 2 based on the performance measure cond{P*wP). 



N 


32 


64 


128 


256 


512 


Choice 1 


1.379 


1.503 


1.621 


1.731 


1.833 


Choice 2 


1.760 


1.887 


2.001 


2.104 


N/A 



3 Digital Shearlet Transform using Compactly Supported Shear- 
lets 

In this section, we will discuss two implementation strategies for computing shearlet coeffi- 
cients associated with a cone-adapted discrete shearlet system now based on compactly sup- 
ported shearlets, as introduced in Chapter Again, one main focus will be on deriving a 
digitization which is faithful to the continuum setting. 

Recall that in the context of wavelet theory, faithful digitization is achieved by the concept 
of multiresolution analysis, where scaling and translation are digitized by discrete operations: 
Downsampling, upsampling and convolution. In the case of directional transforms however, 
three types of operators: Scaling, translation and direction, need to be digitized. In this section, 
we will pay particular attention to deriving a framework in which each of the three operators is 
faithfully interpreted as a digitized operation in digital domain. Both approaches will be based 
on the following digitization strategies: 

• Scaling and translation: A multiresolution analysis associated with anisotropic scaling 
can be applied for each shear parameter k. 
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• Directionality: A faithful digitization of shear operator 52-j/2^ has to be achieved with 
particular care. 

After stating and discussing the two main obstacles we are facing when considering com- 
pactly supported shearlets in Subsection 3.1 we present the digital separable shearlet transform 
(DSST), which is associated with a shearlet system generated by a separable function along- 



side with discussions on its properties, e.g., its redundancy; see Subsection 3.2 Subsection 3.3 



then presents the digital non-separable shearlet transform (DNST), whose shearlet elements are 
generated by non- separable shearlet generator. 



3.1 Problems with Digitization of Compactly Supported Shearlets 

Compactly supported shearlets have several advantages, and we exemplarily mention superior 
spatial localization and simplified boundary adaptation. However, we have to face the following 
two problems: 

(PI) Compactly supported shearlets do not form a tight frame, which prevents utilization of 
the adjoint as inverse transform. 

(P2) There does not exist a natural hierarchical structure, mainly due to the application of a 
shear matrix, which - unlike for the wavelet transform - does not allow a multiresolution 
analysis without destroying a faithful adaption of the continuum setting. 

Let us now comment on these two obstacles, before delving into the details of the imple- 
mentation in Subsection l3.2l 



3.1.1 Tightness 

Let us first comment on the problem of non-tightness. Letting ((T,),e/ denote a frame for L?{M?') 
- for example, a shearlet frame -, each function / e L^(E^) can be reconstructed from its frame 
coefficients ((/, C7i)),e/ by 

f = Y,{f,<y.)s-'i<yi), 

where S — T^ieii'T^d^i the associated frame operator on L^(M?-), see Chapter IT]. However, 
in case that (ct,),^/ does not form a tight frame, it is in general difficult to explicitly compute 
the dual frame elements (cJ,). 

Nevertheless, it is well known that the inverse frame operator ' can be effectively approx- 
imated using iterative schemes such as the Conjugate Gradient method provided that the frame 
{<Ji)iei has 'good' frame bounds in the sense of their ratio being 'close' to 1, see also lf 24ll . 
Therefore, now focussing on the situation of shearlet frames, we may argue that input data / 
can be efficiently reconstructed from its shearlet coefficients, if we use a compactly supported 
shearlet frame with 'good' frame bounds. In fact, the theoretical frame bounds of compactly 
supported shearlet frames have been theoretically estimated as well as numerically computed 
in lfT9l . These results were derived for the class of 2D separable shearlet generators xj/ already 
described in Chapter 1 1 1, which we briefly recall for the convenience of the reader: 

For positive integers K and L fixed, let the ID lowpass filter otq be defined by 

K(^i)|2 = (cos(7r^0rE (^"^^")(sin(7r^i)r, 

for ^1 G M. Further, define the associated bandpass fiher m\ by 

|mi(^i)|2 = |mo(^i + l/2)|2, ^iGE, 
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and the ID scaling function by 

0i(^i) = n'"o(2-^'^i), <^ieK- 

Using the filter mi and scahng function 01, we now define the 2D scaUng function (j) and sepa- 
rable shearlet generator xj/ by 

0(^i,fe) = 0i(^i)0i(fe) and mi (4^0(^1 (^1)^1(2^2). (3.1) 

In |fT9l . it was shown that compactly supported shearlets i/yjt,;n generated by the shearlet gen- 
erator \j/ form a frame for L^(C)^ with appropriately chosen parameters K and L, where 

C = {^eR':\^2/^i\<l, |^i|>l}. 

This construction is directly extended to construct a cone-adapted discrete shearlet frame for 
l2(M2) (cf. also Chapters 0] and H). 

Table [2] provides some numerically estimated frame bounds in L?{C) for certain choice of 
K and L. It shows that indeed the ratio of the frame bounds of this class of compactly supported 
shearlet frames is sufficient small for utilizing an iterative scheme for efficient reconstruction; 
in this sense the frame bounds are 'good'. 

Table 2: Numerically estimated frame bounds for various choices of the parameters K and L. c\ and 
C2 are the sampling constants in the sampling matrix Mc for translation (see Chapter [1]). 



K 


L 




Cl 


B/A 


39 


19 


0.90 


0.15 


4.1084 


39 


19 


0.90 


0.20 


4.1085 


39 


19 


0.90 


0.25 


4.1104 


39 


19 


0.90 


0.30 


4.1328 


39 


19 


0.90 


0.40 


5.2495 



The frequency covering by compactly supported shearlets \j/j,k,m, 

10(^)1'+! L \YiSj:A,,^)\H\^iSlA„^)\\ 

j>OkeKj 

is closely related to the ratio of frames bounds and, in particular, which areas in frequency 
domain cause a larger ratio. This function is illustrated in Fig. [5] which shows that its upper 
and lower bounds are as expected well controlled. 

3.1.2 Hierarchical Structure 

Let us finally comment on the problem to achieve a hierarchical structuring. To allow fast 
implementations, the data structure of the transform is essential. The hierarchical structure of 
the wavelet transform associated with a multiresolution analysis, for instance, enables a fast 
implementation based on filterbanks. In addition, such a hierarchical ordering provides a full 
tree structure across scales, which is of particular importance for various applications such as 
image compression and adaptive PDE schemes. It is in fact mainly due to this property - and 
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(a) Whole frequency plane, (b) Horizontal cone, (c) Vertical cone. 



Figure 5: Frequency covering by shearlets \ \j/j,k^m\'^- (a) Frequency covering of the entire frequency 
plane, (b) Frequency covering of the horizontal cone, (c) Frequency covering of the vertical cone. 

the unified treatment of the continuum and digital setting - that the wavelet transform became 
an extremely successful methodology for many practical applications. 

From a certain viewpoint, shearlets \l/j^k.m can essentially be regarded as wavelets associ- 
ated with an anisotropic scale matrix when the shear parameter k is fixed. This observation 
allows to apply the wavelet transform to compute the shearlet coefficients, once the shear oper- 
ation is computed for each shear parameter k. This approach will be undertaken in the digital 
formulation of the compactly supported shearlet transform, and, in fact, this approach imple- 
ments a hierarchical structure into the shearlet transform. The reader should note that this 
approach does not lead to a completely hierarchical structured shearlet transform - also com- 
pare our discussion at the beginning of this section -, but it will be sufficient for deriving a fast 
implementation while retaining a faithful digitization. 

3.2 Digital Separable Shearlet Transform (DSST) 

We now describe a faithful digitization of the continuum domain shearlet transform based on 
compactly supported shearlets as introduced in [221, which moreover is highly computationally 
efficient. 

3.2.1 Faithful Digitization of the Compactly Supported Shearlet Transform 

We start by discussing those theoretical aspects which allow a faithful digitization of the shear- 
let transform associated with the shearlet system generated by p.l| i. For this, we will only 
consider shearlets \j/j,k.m for the horizontal cone, i.e., belonging to ^'{y/,c). Notice that the 
same procedure can be applied to compute the shearlet coefficients for the vertical cone, i.e., 
those belonging to ^{\jf,c), except for switching the order of variables. 

To construct a separable shearlet generator iff £ L^(M^) and an associated scaling function 
(j) e L^(M^), let (j) e L^(]R) be a compactly supported ID scaling function satisfying 

= ^(«i)v^^i(2-«i -ni) (3.2) 

for some 'appropriately chosen' filter h - we comment on the required condition below. An 
associated compactly supported ID wavelet Xffi e L^(M) can then be defined by 

L ^(ni)\/20i(2xi-ni), (3.3) 

where again g is an 'appropriately chosen' filter. The selected shearlet generator is then defined 
to be 

V/(xi,X2) = ri(xi)0i(x2), (3.4) 
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and the scaling function by 

^{xi,X2) = 0l(xi)0i(x2). 

Let us comment on whether this is indeed a special case of the shearlet generators defined 
in The Fourier transform of y/ defined in ( |3.4| i takes the form 

V/(^i,^2)=mi(^i/2)0i(<^i/2)0i(fe/2), 

where mi is a trigonometric polynomial whose Fourier coefficients are g{ni). We need to 
compare this expression with the Fourier transform of the shearlet generator yf given in 
which is 

V/(^i,^2)=mi(4^i)^i(2^i)0i(fe), 

with ID scaling function defined in p.2| l. We remark that this later scaling function is 
slightly different defined as in This small adaption is for the sake of presenting a simpler 
version of the implementation; essentially the same implementation strategy as the one we will 
describe can be applied to the shearlet generator given in 

The filter coefficients h and g are required to be chosen so that y/ satisfies a certain de- 
cay condition (cf. |19| of Chapter 11]) to guarantee a stable reconstruction from the shearlet 
coefficients. 

For the signal / e L^(M^) to be analyzed, we now assume that, for 7 > fixed, / is of the 
form 

/W= 1^ /y(«)2-'0(2-'xi-«i,2-'x2-«2). (3.5) 

Let us mention that this is a very natural assumption for a digital implementation in the sense 
that the scaling coefficients can be viewed as sample values of / - in fact fj{n) = f{7r^n) with 
appropriately chosen ^. Now aiming towards a faithful digitization of the shearlet coefficients 
(/, y/j^k.m) for /' = 0, . . . ,7 — 1, we first observe that 

(/, Wj.k.,n) = {f{S2-.iP-ki-)). Wj.oA-)), (3.6) 

and, WLOG we will from now on assume that j/2 is integer; otherwise either \j/2] or [y'/2j 
would need to be taken. Our observation ( |3.6[ ) shows us in fact precisely how to digitize the 
shearlet coefficients (/, '//./t.™)- By ^PPlying the discrete separable wavelet transform associ- 
ated with the anisotropic sampling matrix to the sheared version of the data /(5'2-j/2^(-))- 
This however requires - compare the assumed form of / given in p.5| l - that /(52-;/2^(-)) is 
contained in the scaling space 

Vj = {2^</>(2^ • -ni,2-' • -nz) : (ni,n2) & I?}- 

It is easy to see that, for instance, if the shear parameter is non-integer, this is unfor- 

tunately not the case. The true reason for this failure is that the shear matrix S2-j/if. does not 
preserve the regular grid 2^''l? in Vj, i.e.. 

In order to resolve this issue, we consider the new scaling space V'y^_y/2 j defined by 

We remark that the scaling space ^/_^jy2 j is obtained by refining the regular grid 2^^1? along 
the xi-axis by a factor of 2-'/^. With this modification, the new grid 2^^^^I^TL x 2^^TL is now 
invariant under the shear operator 5'2-j/2jl, since with Q ~ diag(2, 1), 

2--^-^'/2z X 2--'Z = 2-^Q-^l^{l})^2-^Q-^l^{Sk{7L^)) 
= S2-j/2i,{2-^'^/^Zx2-^Z). 
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This allows us to rewrite /(52-j/2^-( )) in p.6| l in the following way. 

Lemma 3.1. Retaining the notations and definitions from this subsection, letting f 2'/^ and *i 
denote the ID upsampling operator by a factor of 2^^^ and the ID convolution operator along 
the xi-axis, respectively, and setting hji2{ni) to be the Fourier coefficients of the trigonometric 
polynomial 

%2(^i)= n E/<«i>-'"'''"'^', (3.7) 

we obtain 

f{S2-j/2i^{x)) ^ E fj{Skn)2-'+'/^U2-'+j/\-ni,2'x2-n2), 

neZ^ 

where 

//(n) = ((//)t2^/2*l/^;72)(n)- 
The proof of this lemma requires the following result, which follows from the cascade 
algorithm in the theory of wavelet. 

Proposition 3.1 ([22]). Assume that 0i and Xj/i e L^(M) satisfy equations p.2[ ) and ( |3.3| l re- 
spectively. For positive integers ji < 72, we then have 

2T01 (2/1x1 -m) - E /2,,_;,(^!'i-2>2-^'«i)2T0i (2^2x1-^/1) (3.8) 

diez 

and 

2Tv,i(2>ixi-ni)= gj^^j^{di~V^-^^ni)2'i ^,{V^xi~di), (3.9) 

where hj and gj are the Fourier coefficients of the trigonometric polynomials Hj defined in 
and Gj defined by 

Gji^O = (n L h{n,)e-'^-''"'^^) ( J: gMe-'^'^'-'"'^^) 

for j > fixed. 

Proof of Lemma [377] Equation ( |3.8| l with ji ~ J and j2— J + j/2 implies that 

2-'/20i(2^xi -ni) = ^ V,y2(^i -2^/2«i)2'/2+;74^^(2^+>/2^^ (3 JO) 

Also, since is a 2D separable function of the form 0(xi ,^2) = 0i (^2), we have that 
/W= £ (£ //(«i,«2)2^'"0i(2^^i-«i))2'/2(/)i(2^X2-«2)- 

«2GZ rti^Z 

By ( |3.10| l, we obtain 

/W- 1: /,(«)2^+^'/V(2'e^'/'x-«), 

where Q = diag(2, 1). Using Q-'^^S2-j/2f, — SkQ^I^, this finally implies 

f{S2-,n,{x)) = l^fj{n)2'+^l'H2'&'^S2-,2,{x)-n) 

= £ fj{n)2'+^l',^{Sd2'Q^'^x-S^,n)) 

nGZ2 

= fj{Skn)2'+^l'<^{S,{2-'Q^I'-x-n)). 

nel? 

The lemma is proved. □ □ 
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Figure 6: Illustration of application of the digital shear operator 5^ The dashed lines correspond 
to the refinement of the integer grid. The new sample values lie on the intersections of the sheared 
lines associated with 51/4 with this refined grid. 



The second term to be digitized in p.6[ l is the shearlet i/y itself. A direct corollary from 
Proposition 3.1 is the following result. 

Lemma 3.2. Retaining the notations and definitions from this subsection, we obtain 
Vj,k.m (x) - SJ-j {di - 2'-'nn )hj_j/2 {di - 2'-'l^m2)2'+'l^<^ (2^x - d) . 

As akeady indicated before, we will make use of the discrete separable wavelet transform 
associated with an anisotropic scaling matrix, which, for j\ and y'2 > as well as c € ^(Z^), we 
define by 

^;iJ2('^)("i'"2) = 8jii'ni-2^^ni)hj^_{m2~2->^n2)c{mi,m2), [n\,n2)el?. (3.11) 
Finally, Lemmata 3. 1 and |3.2| yield the following digitizable form of the shearlet coefficients 

{.f,Wj.k.,m)- 

Theorem 3.1 ( 0221 ). Retaining the notations and definitions from this subsection, and letting^ 
2-'/^ be ID downsampling by a factor of 2^1'^ along the horizontal axis, we obtain 

if, Wj.k,m) = "^j-jj-iii ( ) * *i ^72) ^2,/2) ("^)' 

where ^k{n) — {i'{Sk{ ))A{' — n)) for n G 1?', and hji2[n\) — hj/2{—ni). 



3.2.2 Algorithmic Realization 



Computing the shearlet coefficients using Theorem 3.1 now restricts to applying the discrete 
separable wavelet transform ( |3.1 l| l associated with the sampling matrix A^j to the scaling coef- 
ficients 

Before we state the explicit steps necessary to achieve this, let us take a closer look at the scaling 
coefficients S'^'-j/ik^-M' which can be regarded as a new sampling of the data fj on the integer 
grid by the digital shear operator S'^-j/ij^- This procedure is illustrated in Fig. |6jin the case 

2-}/^k = -l/4. 

Let us also mention that the filter coefficients ^^(n) in ( |3.12[ ) can in fact be easily precom- 
puted for each shear parameter k. For a practical implementation, one may sometimes even skip 
this additional convolution step assuming that — X{q.q) ■ 

Concluding, the implementation strategy for the DSST cascades the following steps: 
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• Step 1: For given input data fj, apply the ID upsampling operator by a factor of 2'/^ at 
the finest scale j = J. 

• Step 2: Apply ID convolution to the upsampled input data fj with ID lowpass filter hji2 
at the finest scale j = J. This gives fj. 

• Step 3: Resample fj to obtain fj{Sk{n)) according to the shear sampling matrix at the 
finest scale j = J. Note that this resampling step is straightforward, since the integer grid 
is invariant under the shear matrix S^- 

• Step 4: Apply ID convolution to fj{Sk{n)) with hji2 followed by ID downsampling by 
a factor of 2^/^ at the finest scale j — J. 

• Step 5: Apply the separable wavelet transform Wj_j j_j /2 across scales y = 0,1,.. .,7— 1. 

3.2.3 Digital Realization of Directionality 

Since the digital realization of a shear matrix -Sj-z/a^- by the digital shear operator -Sj-i/z^ 
crucial for deriving a faithful digitization of the continuum domain shearlet transform, we will 
devote this subsection to a closer analysis. 

We start by remarking that in fact in the continuum domain, at least two operators exist 
which naturally provide directionality: Rotation and shearing. Rotation is a very convenient 
tool to provide directionality in the sense that it preserves important geometric information 
such as length, angles, and parallelism. However, this operator does not preserve the integer 
lattice, which causes severe problems for digitization. In contrast to this, a shear matrix does 
not only provide directionality, but also preserves the integer lattice when the shear parameter 
k is integer Thus, it is conceivable to assume that directionality can be naturally discretized by 
using a shear matrix S^- 

To start our analysis of the relation between a shear matrix S2-ji2^- and the associated digital 
shear operator S'^-j^j^^ l^t us consider the following simple example: Set = X{x:xi=o}- Then 
digitize fc to obtain a function fd defined on by setting fd{n) ~ fc{n) for all n £ 1? . For fixed 
shear parameter i G M, apply the shear transform to fc yielding the sheared function (5j ( • ) ) . 
Next, digitize also this function by considering /c('5.s(0)lz2- The functions fd and /c(5i(-))|^2 
are illustrated in Fig. [t] for i = — 1/4. We now focus on the problem that the integer lattice is 




Figure 7: (a) Original image fd{n). (b) Sheared image /c(5'_i/4«). 

not invariant under the shear matrix ^1/4. This prevents the sampling points 81/4(11), n £ I? 
from lying on the integer grid, which causes aliasing of the digitized image fc{S -i/4{-))\j2 as 
illustrated in Fig. |8ja). In order to avoid this aliasing effect, the grid needs to be refined by a 
factor of 4 along the horizontal axis followed by computing sample values on this refined grid. 

More generally, when the shear parameter is given by 5 = — 2^^/^A;, one can essentially 
avoid this directional aliasing effect by refining a grid by a factor of 2-'/^ along the horizontal 
axis followed by computing interpolated sample values on this refined grid. This ensures that 
the resulting grid contains the sampling points ((2^-'/^A;)n2,«2) for any ^2 G ^ and is preserved 
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by the shear matrix S_2-i/ii^- This procedure precisely coincides with the application of the 
digital shear operator S'^-j/?/^^ i-e-> we just described Steps 1-4 from Subsection 3.2.2 in which 
the new scaling coefficients S'^_-p {fj){ri) are computed. 



Let us come back to the exemplary situation of fc — X{x 



-0} 



and 5'_ 



1/4 



we started our 



excursion with and compare fc{S _i/4{-))\^2 with 5 (/j) 1^2 obtained by applying the digital 
shear operator S''_^^^ to /^/. And, in fact, the directional aliasing effect on the digitized image 
fc{S-i/4{n)) in frequency illustrated in Fig. |8ja) is shown to be avoided in Fig. [8](b)-(c) by 
considering J ^4 (/^) I 



Thus apphcation of the digital shear operator S'^_jpi^ allows a faithful 





Figure 8: (a) Aliased image: DFT of fc{S_ 1/4(11)). (b) De- aliased image: S'[^^^{fd){n). (c) De- 
aliased image: DFT of S'^_^^^{fd){n). 

digitization of the shearing operator associated with the shear matrix S2-iiij^. 
3.2.4 Redundancy 

One of the main issues which practical applicability requires is controllable redundancy. To 
quantify the redundancy of the discrete shearlet transform, we assume that the input data / is a 
finite hnear combination of translates of a 2D scaling function (p at scale J as follows; 

n\ =0n2=0 



as it was already the hypothesis in p.5| l. The redundancy - as we view it in our analysis - 
is then given by the number of shearlet elements necessary to represent /. Furthermore, to 
state the result in more generality, we allow an arbitrary sampling matrix Mc — diag(ci ,C2) for 
translation, i.e., consider shearlet elements of the form 



Wj,k,m{-) ^ '2-^W{SkA2j ■ -Mcin). 
We then have the following result. 
Proposition 3.2 ( 11221 '). The redundancy of the DSST is 



1 

C1C2 



Proof. For this, we first consider shearlet elements for the horizontal cone for a fixed scale 
7 € {0, ... ,7— 1}. We observe that there exist 2-'/^+' shearing indices k and 2' • 2'/^ • (ciC2)^' 
translation indices associated with the scaling matrix and the sampling matrix M^, respec- 
tively. Thus, 2^^+' (ciC2)^' shearlet elements from the horizontal cone are required for repre- 
senting /. Due to symmetry reasons, we require the same number of shearlet elements from 
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the vertical cone. Finally, about c^^ translates of the scaling function (j) are necessary at the 
coarsest scale j = 0. 

Summarizing, the total number of necessary shearlet elements across all scales is about 

The redundancy of each shearlet frame can now be computed as the ratio of the number of 
coefficients d„ and this number. Letting 7 — > 0° proves the claim. □ □ 

As an example, choose a translation grid with parameters ci — I and C2 ~ 0.4. Then the 
associated DSST has asymptotic redundancy 10/3. 



3.2.5 Computational Complexity 



A further essential characteristics is the computational complexity (see also Subsection 4.6 1, 
which we now formally compute for the discrete shearlet transform. 

Proposition 3.3 ( 11191 ). The computational complexity of the DSST is 



Proof. Analyzing Steps 1-5 from Subsection 3.2.2 we observe that the most time consuming 



step is the computation of the scaling coefficients in Steps 1 - 4 for the finest scale j — J. This 
step requires ID upsampling by a factor of 2'/^ followed by ID convolution for each direction 
associated with the shear parameter k. Letting L denote the total number of directions at the 
finest scale j — J, and the size of 2D input data, the computational complexity for computing 
the scaling coefficients in Steps 1 - 4 is 0{2^I'^L-N). The complexity of the discrete separable 
wavelet transform associated with for Step 5 requires 0{N) operations, wherefore it is 
negligible. The claim follows from the fact that L = 2(2 •2-'/2 + l). □ □ 

It should be noted that the total computational cost depends on the number L of shear pa- 
rameters at the finest scale j = J, and this total cost grows approximately by a factor of as 
L is increased. It should though be emphasized that L can be chosen in such a way that this 
shearlet transform is favorably comparable to other redundant directional transforms with re- 
spect to running time as well as performance. A reasonable number of di recti ons at the finest 



scale is 6, in which case the constant factor 2'°S2(V2(i/2-i)) jn Proposition 3.3 equals 1. Hence 
in this case the running time of this shearlet transform is only about 6 times slower than the dis- 
crete orthogonal wavelet transform, thereby remains in the range of the running time of other 
directional transforms. 



3.2.6 Inverse DSST 

In Subsection |3.L1| we already discussed that this transform is not an isometry, wherefore the 
adjoint cannot be used as an inverse transform. However, the 'good' ratio of the frame bounds 
in the sense as detailed in Subsection 3. L 1 leads to a fast convergence rate of iterative methods 
such as the conjugate gradient method. Let us mention that using the conjugate gradient method 
basically requires computing the forward DSST and its adjoint, and we refer to li24J and also 
Subsection l2.4.2l for more details. 
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3.3 Digital Non- Separable Shearlet Transform (DNST) 



In this section, we describe an alternative approach to derive a faithful digitization of a discrete 
shearlet transform associated with compactly supported shearlets. This algorithmic realization, 
which was developed in [23 J . resolves the following drawbacks of the DSST: 

• Since this transform is not based on a tight frame, an additional computational effort is 
necessary to approximate the inverse of the shearlet transform by iterative methods. 

• Computing the interpolated sampling values in ( |3.12| l requires additional computational 
costs. 

• This shearlet transform is not shift-variant, even when downsampling associated with 
is omitted. 

We emphasize that although this alternative approach resolves these problems, the algorithm 
DSST provides a much more faithful digitization in the sense that the shearlet coefficients can 
be exactly computed in this framework. 

The main difference between DSST and DNST will be to exploit non-separable shearlet 
generators, which give more flexibility. 



3.3.1 Shearlet Generators 

We start by introducing the non-separable shearlet generators utihzed in DNST. First, define 
the shearlet generator y/™" by 

where the trigonometric polynomial P is a 2D fan filter (c.f. llT2i ). For an illustration of P we 
refer to Fig. |9ta). This in turn defines shearlets V^y™™ generated by non-separable generator 
functions y/" for each scale index j > by setting 

where M^.j is a sampling matrix given by M^.^ = diag(cj , Cj) and Cj and are sampling constants 
for translation. 

One major advantage of these shearlets V'y™,, is the fact that a fan filter enables refinement 
of the directional selectivity in frequency domain at each scale. Fig. [9|a)-(b) show the refined 
essential support of 1//™",, as compared to shearlets y/j,k,m arising from a separable generator as 
in Subsection 13. 2. II 



(b) 




Figure 9: (a) Magnitude response of 2D fan filter, (b) Non-separable shearlet ^"m- ('^) Separable 
shearlet Yj.k,m- 
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3.3.2 Algorithmic Realization 



Next, our aim is to derive a digital formulation of the shearlet coefficients (/, V'y™,,) for ^ 
function / as given in p.5| l. We will only discuss the case of shearlet coefficients associated 
with and S^, the same procedure can be applied for and Si^ except for switching the 
order of variables xi and X2- 



In Subsection 3.2.3 we discretized a sheared function f{S2-j/2j^-) using the digital shear 
operator S'^^_-p^ as defined in p.l2| i. In this implementation, we walk a different path. We digi- 
talize the shearlets yf]Tm(') ~ Vj°0mi^2-J/^k') combining multiresolution analysis and digital 
shear operator S^l^-j/i^. to digitize the wavelet V'^om shear operator ^j-;/^^,, respectively. 

This yields digitized shearlet filters of the form 

Wikin) = S'[-ji2k{pj-jl2 * ^.z) («)' 

where Wj is the 2D separable wavelet filter defined by w/(ni,n2) — gj-j{n\)- hj-j/ii^i) and 
Pj-j/ii.^) the Fourier coefficients of the 2D fan filter given by f (2''^^^i,2''^^''^+'^2)- The 
DNST associated with the non-separable shearlet generators \//"°" is then given by 

DNSTj^,{fj){n) - v7^,)(2-'-^c{«i,2''-^/24«2), for// G 
We remark that the discrete shearlet filters y/^^ are computed by using a similar ideas as in 



Subsection 3.2.1 As before, those filter coefficients can be precomputed to avoid additional 
computational effort. 

Further notice that by setting c\ = 2^^^ and C2 = 2-'/^^"', the DNST simply becomes a 2D 
convolution. Thus, in this case, DNST is shear invariant. 



3.3.3 Inverse DNST 

In case that c\ = 2^^^ and — 2^l^^\ the dual shearlet filters i//^^, can be easily computed by. 



] and we obtain the reconstruction formula 

fj = I,{fj*Wj,k)*vj,k- 

This reconstruction is stable due to the existence of positive upper and lower bounds of Y,j,k W'jki 
which is implied by the frame property of nonseparable shearlets yf"°"„j, see f23l. Thus, no iter- 
ative methods are required for the inverse DNST. The frequency response of a discrete shearlet 



filter Xj/jj^ and its dual xffjj, is illustrated in Fig. 10 We observe that primal and dual shearlet 
filters behave similarly in the sense that both of filters are very well localized in frequency. 



4 Framework for Quantifying Performance 

We next present the framework for quantifying performance of implementations of directional 
transforms, which was originally introduced in |20, 14J. This set of test measures was designed 
to analyze particular well-understood properties of a given algorithm, which in this case are the 
desiderata proposed at the beginning of this chapter. This framework was moreover introduced 
to serve as a tool for tuning the parameters of an algorithm in a rational way and as an objective 
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Figure 10: Magnitude response of shearlet filter y/^^ and its dual filter xj/jj^ 



common ground for comparison of different algorithms. The performance of the three algo- 
rithms FDST, DSST, and DNST will then be tested with respect to those measures. This will 
give us insight into their behavior with respect to the analyzed characteristics, and also allow a 
comparison. However, the test values of these three algorithms will also show the delicateness 
of designing such a testing framework in a fair manner, since due to the complexity of algorith- 
mic realizations it is highly difficile to do each aspect of an algorithm justice. It should though 
be emphasized that - apart from being able to rationally tune parameters - such a framework of 
quantifying performance is essential for an objective judgement of algorithms. The codes of all 
measures are available in lShearLabI 

In the following, S shall denote the transform under consideration, S* its adjoint, and, if 
iterative reconstruction is tested, GaJ shall stand for the solution of the matrix problem AI — J 
using the conjugate gradient method with residual error set to be 10^^. Some measures apply 
specifically to transforms utilizing the pseudo-polar grid, for which purpose we introduce the 
notation P for the pseudo-polar Fourier transform, w shall denote the weighting applied to the 
values on the pseudo-polar grid, and W shall be the windowing with additional 2D iFFT. 



4.1 Algebraic Exactness 

We require the transform to be the precise implementation of a theory for digital data on a 
pseudo-polar grid. In addition, to ensure numerical accuracy, we provide the following measure. 
This measure is designed for transforms utilizing the pseudo-polar grid. 

Measure 1. Generate a sequence of 5 (of course, one can choose any reasonable integer other 
than 5) random images . . J$ on a pseudo-polar grid for N = 512 and R = 8 with standard 
normally distributed entries. Our quality measure will then be the Monte Carlo estimate for the 
operator norm || W*W — Id\\op given by 



Malg ■ 



max - 

1=1... ..5 



\W*WIi-Ii\\2 

Whh ■ 



This measure applies to the FDST - not to the DSST or DNST - , for which we obtain 

Maig = 6.6E-\6. 

This confirms that the windowing in the FDST is indeed up to machine precision a Parseval 



frame, which was already theoretically confirmed by Theorem 2.2 



4.2 Isometry of Pseudo-Polar Transform 

We next test the pseudo-polar transform itself which might be used in the algorithm under 
consideration. For this, we will provide three different measures, each being designed to test a 
different aspect. 
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Measure 2. • Closeness to isometry. Generate a sequence of 5 random images Ii,. . . ,1^ 
of size 512 x 512 with standard uniformly distributed entries. Our quality measure will 
then be the Monte Carlo estimate for the operator norm \\P*wP — Id\\op given by 

\\P*wPIi-Ii\\2 

Misonn = max . 

i=l,-,5 \\Ii\\2 

• Quality of preconditioning. Our quality measure will be the spread of the eigenvalues of 
the Gram operator P*wP given by 

_ An.ax(PVf ) 

• Inveitibility. Our quality measure will be the Monte Carlo estimate for the invertibility 
of the operator ^/wP using conjugate gradient method G (residual error is set to be 
IQ-^, here GaJ means solving matrix problem AI = J using conjugate gradient method) 
given by 

\\G^.py^PIi-Ii\\2 

Misom^ = max T—T . 

/=1,...5 ||/i||2 

This measure applies to the FDST - not to the DSST or DNST - , for which we obtain the 
following numerical results, see Table |3] 

Table 3: The numerical results for the test on isometry of the pseudo-polar transform. 











FDST 


9.3E-4 


1.834 


3.3E-7 



The slight isometry deficiency of Misomi ~ 9.9E-4 mainly results from the isometry defi- 
ciency of the weighting. However, for practical purposes this transform can be still considered 
to be an isometry allowing the utilization of the adjoint as inverse transform. 

4.3 Parseval Frame Property 

We now test the overall frame behavior of the system defined by the transform. These measures 
now apply to more than pseudo-polar based transforms, in particular, to FDST, DSST, and 
DNST 

Measure 3. Generate a sequence of 5 random images Ii, . . . ,Is of size 512 x 512 with standard 
uniformly distributed entries. Our quality measure will then twofold: 

• Adjoint transform. The measure will be the Monte Carlo estimate for the operator norm 
1 1 S*S — Id\\ op given by 

\\S*SIi-Ii\\2 

Mtight, = max . 

'=1 5 ||/;||2 

• Iterative reconstruction. Using conjugate gradient G ^,p, our measure will be given by 

\\G^,pW*SU~h\\2 
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Table 4: The numerical results for the test on Parseval property. 





Mtighti 


Mtight2 


FDST 


9.9E-4 


3.8E-7 


DSST 


1.9920 


1.2E-7 


DNST 


0.1829 


5. BE- 16 (with dual filters) 



The following table, Table |4j presents the performance of FDST and DNST with respect to 
these quantitative measures. 

The transform FDST is nearly tight as indicated by the measures Mtightj = 9.9E-4, i.e., 
the chosen weights force the PPFT to be sufficiently close to an isometry for most practical 
purposes. If a higher accurate reconstruction is required, Mjightj = 3.8E-7 indicates that this can 
be achieved by the CG method. As expected, Mtight, — 1.9920 shows that DSST is not tight. 
Nevertheless, the CG method provides with Mtightj = 1 ^E-? a highly accurate approximation of 
its inverse. DNST is much closer to being tight than DSST (see Mtight, = 0.1829). We remark 
that this transform - as discussed - does not requke the CG method for reconstruction. The 
value 5.8-16 was derived by using the dual shearlet filters, which show superior behavior 



4.4 Space-Frequency-Localization 

The next measure is designed to test the degree to which the analyzing elements, here phrased 
in terms of shearlets but can be extended to other analyzing elements, are space-frequency 
localized. 

Measure 4. Let I be a shearlet in a 512 x 512 image centered at the origin (257,257) with 
slope of scale 4, i.e., CT^ q q + CJ^ q o- ^MaZ/fy measure will be four-fold: 

• Decay in spatial domain. We compute the decay rates di, . . . ^dsi2 along lines parallel to 
the y-axis starting from the line [257, : ] and the decay rates c/512, • • •, 0^1024 with x and y 
interchanged. By decay rate, for instance, for the line [257 : 512, 1], we first compute the 
smallest monotone majorant M{x, \ ), x — 251 , . . . ^512 -note that we could also choose an 
average amplitude here or a different 'envelope ' -for the curve 1 ) |, x = 251 , ...,512. 
Then the decay rate is defined to be the average slope of the line, which is a least square 
fit to the curve log(M(x, 1)), x — 251 , . . . ,512. Basedon these decay rates, we choose our 
measure to be the average of the decay rates 

Mdecay, ^ ^i' 

iVZ.^- 1024 

• Decay in frequency domain. Here we intend to check whether the Fourier transform of I 
is compactly supported and also the decay. For this, let I be the 2D-FFT of I and compute 
the decay rates c/,, / = 1, . . . , 1024 as before. Then we define the following two measures: 

- Compactly supportedness. 

maX|„|.|,,|<3|/(M,v)| 

Msupp = 



max„,v|/(M,v)| 
- Decay rate. 

iUZt ;^l_ 
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• Smoothness in spatial domain. We will measure smoothness by the average of local 
Holder regularity. For each {uq,vq), we compute M{u,v) — \I{u,v) — /(mo,vo)|, < 
max{|M — uq\, \v — Vol} < 4. Then the local Holder regularity o;j,(, is the least square 
fit to the curve log(|M(M, v)|). Then our smoothness measure is given by 

1 ^ 



• Smoothness in frequency domain. We compute the smoothness now for I, the 2D-FFT of 
I to obtain the new CCu^v and define our measure to be 

1 ^ 



Let us now analyze the space-frequency localization of the shearlets utilized in FDST, DSST 
and DNST by these measures. The numerical results are presented in Table|5] 



Table 5: The numerical results for the test on space-frequency localization. 





M decay i 


M.upp 


^ decay 2 


^smooth I 




FDST 


-1.920 


5.5E-5 


-3.257 


1.319 


0.734 


DSST 


— OO 


8.6E-3 


-1.195 


0.012 


0.954 


DNST 


— oo 


2.0E-3 


-0.716 


0.188 


0.949 



The shearlet elements associated with FDST are band-limited and those associated with 
DSST and DNST are compactly supported, which is clearly indicated by the values derived 
for Mdecayi^ M^upp, and Mjecay2- would be expected that Mj^cay^ — ~°° for FDST due to 
the band-limitedness of the associated shearlets. The shearlet elements are however defined 
by their Fourier transform on a pseudo-polar grid, whereas the measure Mdecay, is taken after 
applying the 2D-FFT to the shearlets resulting in data on a cartesian grid, in particular, yielding 
a non-precisely compactly supported function. 

The test values for Mi„„_,o,/,j and Msmoorin show that the associated shearlets are more smooth 
in spatial domain for FDST than for DSST and DNST, with the reversed situation in frequency 
domain. 

4.5 Shear Invariance 

Shearing naturally occurs in digital imaging, and it can - in contrast to rotation - be precisely 
realized in the digital domain. Moreover, for the shearlet transform, shear invariance can be 
proven and the theory implies 

We therefore expect to see this or a to the specific directional transform adapted behavior The 
degree to which this goal is reached is tested by the following measure. 

Measure 5. Let I be an 256 x 256 image with an edge through the origin (129, 129) of slope 
0. Given — 1 < 5 < 1, generates an image 1^ := I{Ss-) and let Sj be the set of all possible scales 
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j such that 2h G Z. Our quality measure will then be the curve 

M shear. j= max r-r , scalej&Sj, 

-2i<k,k+2h<2i \\l\\2 

where Cj± is the shearlet coefficients at scale j and shear k. 
We present our results in Table |6] 

Table 6: The numerical results for the test on shear invariance. 





^ shear, 1 


^shear,2 


Mshear,3 


MghearA 


FDST 


1.6E-5 


1.8E-4 


0.002 


0.003 



This table shows that the FDST is indeed almost shear invariant. A closer inspection shows 
thatMj/,g„r,i ^ndMshear,2 are relatively small compared to the measurements with respect to finer 
scales Msi,ear,3 ^nd MshearA- The reason for this is the aliasing effect which shifts some energy 
to the high frequency part near the boundary away from the edge in the frequency domain. 

We did not test DSST and DNST with respect to this measure, since these transforms show 
a different - not included in this Measure 5 - type of shear invariance behavior 



4.6 Speed 

Speed is one of the most fundamental properties of each algorithm to analyze. Here, we test the 
speed up to a size of = 512 which regard as sufficient to computing the complexity. 

Measure 6. Generate a sequence of 5 random images /,-, / = 5, . . . , 9 of size 2' x 2' with standard 
normally distributed entries. Let Si be the speed of the shearlet transform S applied to li. Our 
hypothesis is that the speed behaves like Sj = c ■ (2^')''; 2^' being the size of the input. Let now 
da be the average slope of the line, which is a least square fit to the curve i log(s, ). Let also 
fi be the 2fft applied to li, i — 5, ... ,9. Our quality measure will then be three-fold: 



Complexity. 



The Constant. 



Mspeed 



da 



' 21og2' 



1^ 



Comparison with 2D-FFT. 



J ,=5 /; 



Table [Tjpresents the results of testing FDST, DSST and DNST with respect to these speed 
measures. 

To interpret these results correctly, we remark that the DNST was tested only with test 
images /, for / = 7, ... ,9, since it can not be implemented for small size images. Interestingly, 
the results also show that the 2D FFT based convolution makes DNST comparable to DSST with 
respect to these speed measures, although it is much more redundant than DSST. Finally, the 
results show that FDST is comparable with both DSST and DNST with respect to complexity 
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Table 7: The numerical results for the test on speed. 





^ speed I 


^speed2 


^speedi 


FDST 


1.156 


9.3E-6 


280.560 


DSST 


0.821 


4.5E-3 


88.700 


DNST 


1.081 


9.9E-8 


40.519 



measure Mspeed, ■ From this, it is conceivable to assume that FDST is highly comparable with 
respect to speed for large scale computations. The larger value Mspeed^ = 280.560 appears due 
to the fact that the FDST employs fractional Fourier transforms on an oversampled pseudo-polar 
grid of size. 

4.7 Geometric Exactness 

One major advantage of directional transforms is their sensitivity with respect to geometric 
features alongside with their ability to sparsely approximate those (cf. Chapter |i4J). This 
measure is designed to analyze this property. 

Measure 7. Let /i, . . . ,/8 be 256 x 256 images of an edge through the origin (129, 129) and of 
slope [—1,-0.5,0,0.5, 1] and the transpose of the middle three, and let ci j be the associated 
shearlet coefficients for image li and scale j. Our quality measure will two-fold: 

• Decay of significant coefficients. Consider the curve 

1 ^ 

- V max |c,- j(of analyzing elements aligned with the line)\, scale j, 
8 ,=1 

let d be the average slope of the line, which is a least square fit to log of this curve, and 
define 

Mgeoi = d. 

• Decay of insignificant coefficients. Consider the curve 

1 ^ 

-^max\cij(of all other analyzing elements)], scale j, 

let d be the average slope of the line, which is a least square fit to log of this curve, and 
define 

^geo2 d. 

Table [8] shows the numerical test results for FDST, DSST, and DNST. 

As expected, the decay rate of the insignificant shearlet coefficients of FDST, i.e., the ones 
not aligned with the line singularity, measured by Mgeo2 ~ -2.032 is much larger than the de- 
cay rate of the significant shearlet coefficients measured by Mgg„, « -1.358. Notice that this 
difference is even more significant in the case of the DSST and DNST. 

4.8 Stability 

To analyze stability of an algorithm, we choose thresholding as the most common impact on a 
sequence of transform coefficients. 
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Table 8 


: The numerical results for the test on 


geometric exactness. 








FDST 


-1.358 


-2.032 


DSST 


-0.002 


-0.030 


DNST 


-0.019 


-0.342 



Measure 8. Let I be the regular sampling of a Gaussian function with mean and variance 
256 on {—128, 127}^ generating an 256 x 256-image. 

• Thresholding 1. Our first quality measure will be the curve 

_ ||G^pW*thresi,p.5/-/H2 

^""""•l.Pl ~ ||/||2 

where thresi_p| discards 100 • (1 — 2^''' ) percent of the coefficients (pi = [2:2: 10]). 

• Thresholding 2. Our second quality measure will be the curve 

_ \\G^.pW'thres2.p,SI~I\\2 

""""•2.P2 ^ ||/|j2 

where thres2,p2 ^^ts all those coefficients to zero with absolute values below the threshold 
m(l — 2^P-) with m being the maximal absolute value of all coefficients. (p2 — [0.001 : 
0.01:0.041]j 

Table|9]shows that even if we discard 100(1 -2-'") - 99.9% of the FDST coefficients, the 
original image is still well approximated by the reconstructed image. Thus the number of the 
significant coefficients is relatively small compared to the total number of shearlet coefficients. 
From Table [TO] we note that knowledge of the shearlet coefficients with absolute value greater 
than m(l — 1/2"-'"" )(^ 0.1% of coefficients) is sufficient for precise reconstruction. 

DNST shows a similar behavior with worse values for relatively large pi. It should be 
however emphasized that firstly, the redundancy of DNST used in this test is 25 and this is 
lower than the redundancy of FDST, which is about 71. This effect can be more strongly 
seen by the test results of DSST whose redundancy with 4 even much smaller Secondly, a 
significant part of the low frequency coefficients in both DSST and DNST will be removed by a 
relatively large threshold, since the ratio between the number of the low frequency coefficients 
and the total number of coefficients is much higher than FDST. This prohibits a similarly good 
reconstruction of a Gaussian function. 

This test in particular shows the delicateness of comparing different algorithms by merely 
looking at the test values without a rational interpretation; in this case, without considering the 
redundancy and the ratio between the number of the low frequency coefficients and the total 
number of coefficients. 
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PI 2 A 
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6 8 10 

2.5E-05 0.001 0.007 

0.02961 0.0296 0.0331 

0.00391 0.0124 0.0396 

3.4E-05 5.5E-05 5.5E-05 
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P2 


0.001 
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0.021 


0.031 


0.041 


FDST 


0.005 


0.039 


0.078 


0.113 


0.154 


DSST 


0.030 


0.036 


0.046 


0.056 


0.072 


DNST 


0.002 


0.018 
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0.055 


0.076 


DNST+DWT 0.001 


0.013 


0.020 


0.024 


0.031 
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